home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / DESCRIPT.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  8KB  |  314 lines
  1. /* 
  2.  * descript.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * DESCRIPT
  11.  *
  12.  *  evaluation of Fourier descriptors for a plane curve
  13.  *   encoded in a set {delta_phi(k), l(k)};
  14.  *
  15.  *  references:
  16.  *  C. T. Zahn, 
  17.  *  Proc. Joint Intern. Conf. on Artif. Intelligence, May 1969, pp. 621 - 628;
  18.  *  SLAC Report 70, Stanford Linear Accelerator Center, 1968;
  19.  *  C. T. Zahn & R. Z. Roskies, IEEE Trans. Comp., Vol C-21, 269 - 281 (1972);
  20.  *
  21.  *  -->the truncated Fourier expansion derived by Z & R is computed 
  22.  *     in form of inner products using array processor functions
  23.  *
  24.  *  -->test shapes (see Z & R): 
  25.  *           four (xdrawfour.c), inverted L (fdzahnl.c)
  26.  *
  27.  * ---------------------------------------------------------------------------
  28.  *
  29.  */
  30. #include <stdio.h>
  31. #include <stdlib.h>
  32. #include <math.h>
  33. #include "ip.h"
  34.  
  35. #define    EPS        0.00000001
  36. #define    SQ2        1.414213562
  37.  
  38. #define    ON        1
  39. #define    OFF        0
  40. #define    FD_DEBUG    ON
  41.  
  42.  
  43. /*
  44.  * descriptors()
  45.  *   DESCRIPTION:
  46.  *     evaluate ZAHN-ROSKIES Fourier descriptors
  47.  *   ARGUMENTS:
  48.  *     dphik: delta phik for polygon (float *)
  49.  *     dlk: delta lk for polygon (float *)
  50.  *     mcp: length of dphik and dlk (long)
  51.  *     a_n: Fourier coefficients array (float *)
  52.  *     b_n: Fourier coefficients array (float *)
  53.  *     n_order: order of Fourier descriptors to compute (long)
  54.  *   RETURN VALUE:
  55.  *     none
  56.  */
  57.  
  58. void
  59. descriptors (float *dphik, float *dlk, long mcp, float *a_n, float *b_n, long n_order)
  60. {
  61.   int k, n_ord;
  62.   int fund_flag = OFF;
  63.  
  64.   float real_length, length;
  65.   float dc_term, a, b, part_sum;
  66.   float *angle;
  67.   float *sne, *csn;
  68.   float *mag, *phase;
  69.   int n;
  70.  
  71.  
  72. /*  
  73.  *  compute (in place) running sums of the elements in dlk 
  74.  *  and the corresponding phase angles required in the evaluation
  75.  *  of the fundamental
  76.  */
  77.   if ((sne = (float *) calloc ((int) mcp, sizeof (double))) == NULL) {
  78.     printf ("descriptors: memory allocation error for sne\n");
  79.     exit (1);
  80.   }
  81.   if ((csn = (float *) calloc ((int) mcp, sizeof (double))) == NULL) {
  82.     printf ("descriptors: memory allocation error for csn\n");
  83.     exit (1);
  84.   }
  85.   if ((angle = (float *) calloc ((int) mcp, sizeof (double))) == NULL) {
  86.     printf ("descriptors: memory allocation error for angle\n");
  87.     exit (1);
  88.   }
  89.  
  90.  
  91. /*
  92.  * evaluate contour length, based on curvature points
  93.  */
  94.   part_sum = *(dlk + 0);
  95.   for (k = 1; k < mcp; k++) {
  96.     part_sum += *(dlk + k);
  97.     *(dlk + k) = part_sum;      /* array of partial sums */
  98.   }
  99.   length = *(dlk + mcp - 1);
  100.  
  101.   for (k = 0; k < mcp; k++)
  102.     *(angle + k) = (float) (2 * PI * (*(dlk + k)) / length);
  103.  
  104.   real_length = (float) (0.5 * length);
  105.   printf ("\n...length = %f\n", real_length);
  106.  
  107. /*
  108.  *  dphik contains elements of delta_phik multiplied by PI/4;
  109.  *  to speed up computation, factor out PI 
  110.  */
  111.   for (k = 0; k < mcp; k++)
  112.     *(dphik + k) *= (float) 0.25;
  113.  
  114. #ifdef FD_DEBUG
  115.   printf ("...elements of dphik and dlk:\n");
  116.   for (k = 0; k < mcp; k++)
  117.     printf ("%7.4f  %7.4f\n", *(dphik + k), *(dlk + k));
  118. #endif
  119.  
  120. /*
  121.  *  compute DC term
  122.  */
  123.   vdot (dphik, dlk, &dc_term, mcp);
  124.   dc_term = (float) (-PI * (1.0 + dc_term / length));
  125.   printf ("...and we have a result for the DC term: %f   \n", dc_term);
  126.  
  127. /*
  128.  * compute increasing orders of Fourier coefficients up to n
  129.  */
  130.   n = 0;
  131.   while (n < n_order) {
  132.     for (k = 0; k < mcp; k++) {
  133.       if (fabs (*(sne + k) = (float) sin (*(angle + k) * (n + 1))) < EPS)
  134.         *(sne + k) = (float) 0.0;
  135.       if (fabs (*(csn + k) = (float) cos (*(angle + k) * (n + 1))) < EPS)
  136.         *(csn + k) = (float) 0.0;
  137.     }
  138.     vdot (dphik, sne, &a, mcp);
  139.     vdot (dphik, csn, &b, mcp);
  140.     *(b_n + n) = (b / (float) (n + 1));
  141.     *(a_n + n) = -(a / (float) (n + 1));
  142.  
  143.     n++;
  144.   }
  145.   printf ("...descriptors up to %dth order are:\n", n_order);
  146.   for (n = 0; n < n_order; n++)
  147.     printf ("%7.3f   %7.3f\n", *(a_n + n), *(b_n + n));
  148.  
  149.   mag = (float *) calloc ((int) n_order, sizeof (float));
  150.   phase = (float *) calloc ((int) n_order, sizeof (float));
  151. /*
  152.  *  calculate xy to polar coordinate transformation
  153.  */
  154.   vrtop (a_n, b_n, mag, phase, n_order);
  155.   printf ("...polar descriptors up to %ld-th order are:\n", n_order);
  156.   for (n_ord = 0; n_ord < (int) n_order; n_ord++) {
  157.     if (*(phase + n_ord) < 0.0)
  158.       *(phase + n_ord) += (float) (2.0 * PI);
  159.     printf ("%7.3f   %7.3f\n", *(mag + n_ord), *(phase + n_ord));
  160.   }
  161.  
  162.   free (angle);
  163.   free (sne);
  164.   free (csn);
  165.   free (mag);
  166.   free (phase);
  167. }
  168.  
  169. /*
  170.  * vdot()
  171.  *   DESCRIPTION:
  172.  *     take the dot product of 2 vectors
  173.  *     and compute the sum
  174.  *   ARGUMENTS:
  175.  *     v1: first vector (float *)
  176.  *     v2: second vector (float *)
  177.  *     vsum: vector sum (float *)
  178.  *     vlen: length of vectors (long)
  179.  *   RETURN VALUE:
  180.  *     none
  181.  */
  182.  
  183. void
  184. vdot (float *v1, float *v2, float *vsum, long vlen)
  185. {
  186.   float rsum = (float) 0.0;
  187.   int i;
  188.  
  189.   for (i = 0; i < vlen; i++)
  190.     rsum = rsum + *(v1 + i) * (*(v2 + i));
  191.   *vsum = rsum;
  192. }
  193.  
  194. /*
  195.  * vcos()
  196.  *   DESCRIPTION:
  197.  *     takes the cosine of a vector
  198.  *   ARGUMENTS:
  199.  *     source: input vector (float *) NOTE: radians!!
  200.  *     dest: output vector (float *)
  201.  *     vlen: vector length (long)
  202.  *   RETURN VALUE:
  203.  *     none
  204.  */
  205.  
  206. void
  207. vcos (float *source, float *dest, long vlen)
  208. {
  209.   int i;
  210.  
  211.   for (i = 0; i < vlen; i++)
  212.     *(dest + i) = (float) cos ((double) *(source + i));
  213. }
  214.  
  215. /*
  216.  * vsin()
  217.  *   DESCRIPTION:
  218.  *     takes the sine of a vector
  219.  *   ARGUMENTS:
  220.  *     source: input vector (float *) NOTE: radians!!
  221.  *     dest: output vector (float *)
  222.  *     vlen: vector length (long)
  223.  *   RETURN VALUE:
  224.  *     none
  225.  */
  226.  
  227. void
  228. vsin (float *source, float *dest, long vlen)
  229. {
  230.   int i;
  231.  
  232.   for (i = 0; i < vlen; i++)
  233.     *(dest + i) = (float) sin ((double) *(source + i));
  234. }
  235.  
  236. /*
  237.  * vrtop()
  238.  *   DESCRIPTION:
  239.  *     converts rectangular cartesian
  240.  *     coordinates to polar
  241.  *   ARGUMENTS:
  242.  *     x: input x vector (float *)
  243.  *     y: input y vector (float *)
  244.  *     mod: modulus output vector (float *)
  245.  *     arg: atan2 output vector (float *)
  246.  *     vlen: vector lengths (long)
  247.  *   RETURN VALUE:
  248.  *     none
  249.  */
  250.  
  251. void
  252. vrtop (float *x, float *y, float *mag, float *phase, long vlen)
  253. {
  254.   int i;
  255.  
  256.   for (i = 0; i < vlen; i++) {
  257.     *(mag + i) = (float) sqrt ((*(x + i)) * (*(x + i)) + (*(y + i)) * (*(y + i)));
  258.     *(phase + i) = (float) atan2 (*(y + i), *(x + i));
  259.     if (fabs (*(phase + i)) < EPS)
  260.       *(phase + i) = (float) 0.0;
  261.     else if (*(phase + i) < (float) 0.0)
  262.       *(phase + i) += (float) (2.0 * PI);
  263.   }
  264. }
  265.  
  266. /*
  267.  * msdescriptors()
  268.  *   DESCRIPTION:
  269.  *     wrapper for call to descriptors()
  270.  *   ARGUMENTS:
  271.  *     delta_phik: delta phik for polygon (float *)
  272.  *     delta_lk: delta lk for polygon (float *)
  273.  *     mcp: length of dphik and dlk (long)
  274.  *     a_n: Fourier coefficients array (float *)
  275.  *     b_n: Fourier coefficients array (float *)
  276.  *     n_order: order of Fourier descriptors to compute (long)
  277.  *   RETURN VALUE:
  278.  *     none
  279.  */
  280.  
  281. void
  282. msdescriptors (float *delta_phik, float *delta_lk, long mcp, float *a_n, float *b_n, long n_order)
  283. {
  284.   int k;
  285.  
  286. #ifdef FD_DEBUG
  287.   printf ("\nevaluation of Fourier descriptors (Zahn & Roskies)\n");
  288.   printf ("--------------------------------------------------\n");
  289.  
  290.   printf ("...given: two vectors, delta_phik[] and delta_lk[]\n");
  291.  
  292.   printf (" delta_phik:\n");
  293.   for (k = 0; k < (int) mcp; k++) {
  294.     printf ("%7.2f  ", *(delta_phik + k));
  295.     if ((k + 1) % 8 == 0)
  296.       printf ("\n");
  297.   }
  298.   printf ("\n   delta_lk:\n");
  299.   for (k = 0; k < (int) mcp; k++) {
  300.     printf ("%7.2f  ", *(delta_lk + k));
  301.     if ((k + 1) % 8 == 0)
  302.       printf ("\n");
  303.   }
  304. #endif
  305.  
  306.   if (n_order > mcp) {
  307.     printf ("...polygon only has %ld vertices;");
  308.     printf (" cannot yield %ld harmonics!\n", mcp, n_order);
  309.     exit (1);
  310.   }
  311.   descriptors (delta_phik, delta_lk, mcp, a_n, b_n, n_order);
  312.  
  313. }
  314.